source("read.data.R")
source("util.R")
source("analy.model.R")
source("analy.model.selection.R")
source("analy.QRE.general.R")
source("analy.plot.R")

library(lme4)
library(lmerTest)

# e0<-par.setting[1] 
# ea<-par.setting[2]
# inv.cost<-par.setting[3]
# p0<-par.setting[4]
# p1<-par.setting[5]
# penalty<-par.setting[6]
# 
# gamma.pay<-par.behavior[1]
# gamma.r<-par.behavior[2]
# gamma.a<-par.behavior[3]
# gamma.i<-par.behavior[4]

time.start<-proc.time()

library(parallel)
cluster<-makeCluster(6)
#clusterExport(cl=cluster,lsf.str(),envir=.GlobalEnv)
#clusterExport(cl=cluster,ls(),envir=environment())

#treat.list<-c("Base","Penalty","Chat","Norm","Should_invest","Should_not_pay","73%_Invest","62%_Not_Pay","38%_Pay","Chat_Penalty")
#treat.list<-c("Base","Should_invest","Should_not_pay","73%_Invest","62%_Not_Pay","Penalty","Chat","Chat_Penalty")
#treat.list<-c("Base","Should_invest","Should_not_pay")
treat.list<-c("Base")
# treat.list<-c("Should_invest","Should_not_pay")
#treat.list<-c("62%_Not_Pay")
#ff.model.list<-c("analy.model.bias.att.prob.r","analy.model.bias.inv.cost.R")
#model.list<-c("att-prob-bias","inv-cost-bias")
ff.model.list<-c("analy.model.R")
model.list<-c("final model")

clusterExport(cl=cluster,ls(),envir=.GlobalEnv)

# cn1<-c("gamma pay","gamma ransom","gamma attack","gamma invest","prob att bias","pay bias","prev att bias")
# cn2<-c("gamma pay","gamma ransom","gamma attack","gamma invest","inv cost bias","chat pay impact","pay social norm","inv social norm","chat invest impact")
cn2<-c("gamma pay","gamma ransom","gamma attack","gamma invest","follow me inv","pay social norm","pay fairness","inv social norm","follow you inv")
# cn2<-c("gamma pay","gamma ransom","gamma attack","gamma invest","N/A","pay social norm","pay fairness","inv social norm","loss aversion")

#cn.list<-list(cn1,cn2)
cn.list<-list(cn2)

# data<-read.file.index(1)
# Treatment<-get.data(4,"Treatment")
# data<-cbind.data.frame(data,Treatment)
# attack.decision<-get.data(4,"attack.decision")
# data<-cbind.data.frame(data,attack.decision)

data<-read.file.index(1)
data<-subset(data,data[,"Period"]>1)

#a<-subset(data,data[,"Treatment"]=="Base")
#a<-subset(a,a[,"Group.ID"]==1)

result.summary<-NULL
for(treat in treat.list) for(model.index in 1:length(model.list)) {

if(treat!="Pooled") data.base<-subset(data,data[,"Treatment"]==treat) 
else {
#  data.base<-data
  data.base<-subset(data,data[,"Treatment"]!="Norm")
}
source(ff.model.list[model.index])

#llk.test1<-llk(c(0.032107282,0.131139841,0.049238119,0.128641201,-0.489061094,-11.54054904,12.11456325),
#               par.restrict=c(1,1,1,1,1,1,1,0,0),par.restrict.value=rep(0,9),data=data.base)
#llk.test2<-llk(c(0.1,0.1,0.1,0.0,-0.9,-10),par.restrict=c(1,1,1,1,1,1),par.restrict.value=c(0,0,0,0),data=data.base)

#llk.test1<-llk(c(0.1,0.1,0.1,0.1,-2,0,0,5,5),par.restrict=rep(1,9),par.restrict.value=rep(0,9),data=data.base)
#cat(treat," llk=",llk.test1,"\n")

# browser()

ms<-model.selection(data.base,start=c(0.1,0.1,0.1,0.01,0,0,0,0,0),mode="static") 
#ms<-model.selection(data.base,start=c(0.0246,0.0597,0.0742,0.0448,-0.47,3,12.29,46.54,3),mode="static") 
# ms<-model.selection(data.base,start=c(0.0361,0.102,0.0584,0.01,-0.1,0,35.7,31.0,0),mode="static") 
#ms<-model.selection(data.base,start=c(0.0246,0.0597,0.0742,0.0448,-0.47,3,12.29,46.54,3),mode="static") 
# ms<-model.selection(data.base,start=c(0.05,0.05,0.05,0.1,-0.4,0,6,15,0),mode="static") 

#browser()

cn<-cn.list[model.index]

model.summary<-ms.summary(ms,coef.name=cn)

par.setting<-c(data.base[1,"Data1"],data.base[1,"Outside"],data.base[1,"Investment1"],0.8,0.5,data.base[1,"Penalty"],0,0,0)
par.b<-convert.par.stl(ms[[2]]$par,model.spec=ms[[1]])
par.b<-as.numeric(par.b)
par.b<-ifelse(is.na(par.b),ms[[7]],par.b)
# 
# r<-pmax(1,data.base[,"ransom"])
# inv<-ifelse(data.base[,"attack.decision"]==1,data.base[,"i_1"],data.base[,"i_2"])
# pay<-ifelse(data.base[,"attack.decision"]==1,data.base[,"p_1"],data.base[,"p_2"])
# 
# pp1<-prob.pay(c(1:100),1,par.setting=par.setting,par.behavior=par.b)
# pp0<-prob.pay(c(1:100),0,par.setting=par.setting,par.behavior=par.b)
# #pp.model<-mean(na.omit(pp[r]))
# #pp.fit<-mean(na.omit(c(data.base[,"p_1"],data.base[,"p_2"])))
# 
# data.pay<-cbind(r,inv,pay)
# data.pay<-subset(data.pay,(!is.na(data.pay[,"r"]))&(!is.na(data.pay[,"pay"])))
# 
# pp.data<-c(mean(data.pay[data.pay[,"inv"]==0,"pay"]),mean(data.pay[data.pay[,"inv"]==1,"pay"]))
# pp.model<-c(mean(pp0[data.pay[,"r"]][data.pay[,"inv"]==0]),mean(pp1[data.pay[,"r"]][data.pay[,"inv"]==1]))
# 
# pp.summary<-cbind(c(0,1),pp.model,pp.data)
# dimnames(pp.summary)[[2]]<-c("inv","model","data")
# 
# 
# qrr1<-qre.r(1,par.setting=par.setting,par.behavior=par.b)
# qrr0<-qre.r(0,par.setting=par.setting,par.behavior=par.b)
# 
# qrr<-cbind(qrr1[,1],rep(0,nrow(qrr1)),rep(0,nrow(qrr1)))
# 
# for(i in 1:nrow(data.base)) {
#   if(data.base[i,"attack.decision"]>0) {
#     r<-max(1,data.base[i,"ransom"])
#     
#     old.payment.outcome<-ifelse(is.na(data.base[i,"old_payment_outcome"]),0,data.base[i,"old_payment_outcome"])
#     par.setting[7]<-old.payment.outcome
#     
#     ii<-ifelse(data.base[i,"attack.decision"]==1,data.base[i,"i_1"],data.base[i,"i_2"])
#     qrr.i<-qre.r(ii,par.setting=par.setting,par.behavior=par.b)
#     
#     qrr[r,3]<-qrr[r,3]+1
#     if(r<1) browser()
#     if(r>100) browser()
#     
#     qrr[,2]<-qrr[,2]+qrr.i[,2]
#   }
# }
# 
# qrr[,2]<-qrr[,2]/sum(qrr[,2])
# qrr[,3]<-qrr[,3]/sum(qrr[,3])
# 
# dimnames(qrr)[[2]]<-c("r","model","data")
# 
# qrr.plot<-plot.smooth(qrr,n=10)
# 
# ff<-paste("output/figure.",format(Sys.time(), "%Y_%m_%d_%H_%M_%S"),".",treat,".",model.list[model.index],".jpg",sep="")
# 
# m.lab<-paste(treat," ",model.list[model.index],sep="")
# plot.multiple(qrr.plot[,1],qrr.plot[,c(2,3)],c("model","data"),file=ff,main=m.lab,xlab="ransom",ylab="QRE prob")
# 
att.summary<-cbind(c(0,1,2),c(0,0,0),c(0,0,0))
for(i in 1:nrow(data.base)) {
   inv.list<-c(data.base[i,"i_1"],data.base[i,"i_2"])
   qra<-qre.att(inv.list,par.setting=par.setting,par.behavior=par.b)
   att.summary[,2]<-att.summary[,2]+qra[,2]
   att.summary[(att.summary[,1]==data.base[i,"attack.decision"]),3]<-att.summary[(att.summary[,1]==data.base[i,"attack.decision"]),3]+1
}
 
att.summary[,2]<-att.summary[,2]/sum(att.summary[,2])
att.summary[,3]<-att.summary[,3]/sum(att.summary[,3])
 
dimnames(att.summary)[[2]]<-c("att","model","data")
# 
# 
u.m1<-matrix(nrow=2,ncol=2)
u.m2<-matrix(nrow=2,ncol=2)
 
for(i in 0:1) for(j in 0:1) {
   u.m1[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
   u.m2[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
}
 
dimnames(u.m1)<-list(c("0","1"),c("0","1"))
dimnames(u.m2)<-list(c("0","1"),c("0","1"))
 
#qre.i<-find.qre(c(par.b[4],par.b[4]),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
qre.i<-find.qre(c(0.01,0.01),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")

qi1<-qre.i[[1]]
qi2<-qre.i[[2]]
inv.summary<-c("no-no","no-yes","yes-no","yes-yes")
inv.summary<-cbind(inv.summary,c(qi1[1,2]*qi2[1,2],qi1[1,2]*qi2[2,2],qi1[2,2]*qi2[1,2],qi1[2,2]*qi2[2,2]))
 
i1<-data.base[,"i_1"]
i2<-data.base[,"i_2"]
 
i.count<-c(sum((i1==0)&(i2==0)),sum((i1==0)&(i2==1)),sum((i1==1)&(i2==0)),sum((i1==1)&(i2==1)))
i.freq<-i.count/sum(i.count)
 
inv.summary<-cbind(inv.summary,i.freq)
dimnames(inv.summary)[[2]]<-c("inv","model","data")

ans.i1<-glmer(i_1~old_i_1+old_i_2+(1|Group.ID),data.base,family=binomial)
ans.i2<-glmer(i_2~old_i_1+old_i_2+(1|Group.ID),data.base,family=binomial)

coln<-c("treat/model",dimnames(model.summary)[[2]])
tm.col<-c(treat,model.list[model.index],rep(" ",nrow(model.summary)-2))
model.summary<-cbind.data.frame(tm.col,model.summary)

#pp.summary<-rbind(c("pay prob","pay prob","pay prob"),pp.summary,matrix(rep(" ",3*(nrow(model.summary)-1-nrow(pp.summary))),ncol=3))
#att.summary<-rbind(c("att prob","att prob","att prob"),att.summary,matrix(rep(" ",3*(nrow(model.summary)-1-nrow(att.summary))),ncol=3))
#inv.summary<-rbind(c("inv prob","inv prob","inv prob"),inv.summary,matrix(rep(" ",3*(nrow(model.summary)-1-nrow(inv.summary))),ncol=3))
 
# model.summary<-cbind.data.frame(model.summary,pp.summary,att.summary,inv.summary)
result.summary<-rbind(result.summary,model.summary)
}

stopCluster(cluster)

time.end<-proc.time()
tt<-time.end[3]-time.start[3]
tt<-tt/60
cat("time elapsed = ",tt," minutes\n")
save.image()